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1  Background 

This  final  report  details  the  progress  made  under  the  Aero- Optics  Code  Development  Program  managed  by 
the  Computational  Sciences  Branch,  Air  Vehicles  Directorate,  Air  Force  Research  Laboratory.  This  program 
was  tasked  with  developing  an  in-house  capability  to  perform  computational  analysis  of  the  aero-optical 
environment  associated  with  configuration  of  interest  to  the  Air  Vehicles  Directorate.  The  progress  in  the 
first  two  years  focused  on  two  main  areas:  1)  identify  and  obtain  high  quality  experimental  data  for  various 
aero-optic  configurations  to  serve  as  possible  validation  data  for  future  code  development,  and  2)  evaluate 
the  unstructured  flow  solver  AVUS  for  use  on  canonical  cases  representing  the  types  of  problems  expected  to 
be  encountered  in  the  aerodynamic  analysis  of  future  aero-optics  configurations  and  improve  the  predictive 
capabilities  of  the  code  if  feasible.  The  results  of  the  effort  in  these  two  areas  may  be  found  in  Reference  [1]. 
The  third  area  addressed  during  this  effort  involved  the  development  of  specialized  versions  of  the  CFD 
solvers  OVERFLOW  and  AVUS  to  study  aerodynamically-induced  aberrations  to  optical  wave  forms  as 
they  are  propagated  through  an  unsteady,  compressible  and  turbulent  flow  field.  The  initial  work  in  this 
area  may  be  found  in  Reference  [2],  where  an  aero-optics  capability  was  added  to  the  OVERFLOW  solver. 
The  current  work  details  the  follow-on  effort  to  improve  the  implementation  of  the  aero-optics  solver  in 
OVERFLOW  based  on  lessons  learned,  and  to  extend  the  approach  to  the  unstructured  solver  AVUS. 
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2  Introduction 


The  propagation  of  electromagnetic  beams  through  turbulent  flows  has  been  an  important  topic  for  many 
years,  with  applications  ranging  from  missile  defense,  directed  energy,  to  target  designation  and  tracking. 
An  important  aspect  of  these  applications  is  determining  the  effective  beam-on-target  characteristics  after 
the  beam  has  propagated  through  both  the  near-field  turbulent  flow  field  of  the  vehicle  and  the  far-held 
turbulent  atmosphere.  Near-held  propagation  (i.e.  aero-optics)  maintains  some  similarities  to  the  far- 
held  (i.e.  atmospheric)  propagation,  but  due  to  the  interactions  between  turbulence  length  scales,  beam 
wavelengths,  apertures  and  distances,  the  two  often  require  different  approaches  [3,  4,  5,  6].  As  shown  in 
Sutton  [4],  the  techniques  needed  to  evaluate  the  propagation  of  a  beam  from  the  vehicle  to  the  target,  or 
to  evaluate  the  image  imposed  on  a  sensor  from  the  targets  signature  are  similar,  although  the  engineering 
challenges  from  the  two  are  quite  different  (e.g..;  adaptive  optics  versus  remote  sensing). 

Many  noteworthy  efforts  have  been  made  applying  Computational  Fluid  Dynamics  models  for  aero- 
optics  predictions.  For  example,  Tsai  [7]  in  1990  performed  Direct  Numerical  Simulations  (DNS)  of  a 
two-dimensional,  planar  mixing  layer,  and  evaluated  phase  differences  via  Optical  Path  Difference  (OPD) 
integration.  Large  Eddy  Simulations  (LES)  have  been  made  by  many,  including  Jones  [8],  Sinha  [9]  Mani  [10, 
11]  and  recently  by  Visbal  [12],  and  others.  In  these  calculations,  OPD  variations,  or  equivalently,  phase 
differences,  have  been  used  to  characterize  the  beam  distortion  caused  by  the  unsteady  turbulent  flow  field.  A 
Partially- Averaged  Navier-Stokes  (PANS)  approach  has  been  applied  by  Ceniceros  [13,  14]  to  a  turret  model, 
where  significant  correlation  to  experimental  statistics  was  shown,  by  comparing  the  unsteady  variations  of 
OPD  about  a  mean.  An  approach  which  solves  the  Reynolds- Averaged  Navier-Stokes  (RANS)  equations 
coupled  with  a  compressibility  corrected,  two-equation  turbulence  model  and  a  transport  equation  for  the 
density  perturbation  variances  have  been  made  in  Pond  [15],  where,  again,  OPD  is  used  to  characterize  the 
aero-optical  distortion.  The  approach  of  solving  a  density  perturbation  variance  transport  equation  was  also 
used  previous  to  this  work  by  Smith  [16]  within  a  parabolized/thin-layer  solver,  valid  for  shear  layers. 

The  assumptions  often  made  to  apply  a  purely  geometric  optics  approach,  such  as  those  noted  above, 
may  not  be  met  for  all  flow  and  beam  propagation  conditions.  As  noted  in  Truman  [17],  diffraction  may  be 
neglected  when  XL /l2  «  1  (i.e.:  the  Fresnel  condition),  where  the  beam  wavelength  is  A,  the  propagation 
length  is  L  and  the  turbulence  length  scale  is  l.  The  geometric  optics  approach  is  based  upon  the  reduction 
of  Maxwells  equations  to  a  parabolic  (paraxial)  form,  and  then  further  reducing  this  form  to  the  solution 
of  a  phase-shift  operator  [17].  This  assumes  that  the  phase  shift  has  no  impact  upon  the  beam  amplitude, 
although  the  interaction  between  the  phase  and  the  beam  spreading  may  not  be  negligible  for  all  beam 
forms,  propagation  distances  and  index  of  refraction  profiles.  Since  many  applications  need  to  predict 
effects  such  as  scintillation,  the  application  of  a  purely  geometric  approach  might  not  be  broadly  applicable, 
and  care  needs  to  be  taken  in  determining  when  to  apply  it.  For  shorter  wavelengths,  the  phase  shifting 
is  more  effective  (due  to  the  increased  indices  of  refraction),  and  the  scale  difference  between  turbulence 
and  beams  are  closer,  so  their  interaction  with  the  beam  diffraction  is  more  pronounced.  Assumptions 
regarding  compressibility,  turbulence  equilibrium,  scales  and  isotropy  are  also  often  made  [3,  4,  5,  6]  that 
may  not  be  applicable  to  the  flow  in  the  vicinity  of  an  aerospace  vehicle.  As  noted  in  Siegenthaler  [5]  and 
in  Gordeyev  [18,  19,  20],  Visbal  [12]  and  others,  the  complex  flow  patterns  produced  by  protuberances,  such 
as  turrets,  induce  anisotropic  turbulence  fields,  possibly  not  in  turbulent  equilibrium,  with  a  wide  range  of 
length  scales.  With  an  increase  in  vehicle  speeds,  the  assumptions  often  made  regarding  incompressibility 
are  also  suspect. 

In  Reference  [1,  21],  a  specialized  version  of  the  OVERFLOW  2. lx  overset  CFD  solver  was  developed 
to  study  aerodynamically-induced  aberrations  to  optical  wave  forms  as  they  are  propagated  through  an 
unsteady,  compressible  and  turbulent  flow  field.  A  reduced  form  of  Maxwells  equations,  the  Paraxial  Beam 
equation,  was  solved  concurrently  with  the  Navier-Stokes  equations  to  propagate  a  given  wave  form  through 
the  CFD  generated  density  field,  produced  by  Large  Eddy  Simulation  (LES)  or  Detached  Eddy  Simulation 
(DES).  The  integration  with  OVERFLOW  was  made  through  a  Fortran90  module  containing  procedures 
invoked  by  the  master  OVERFLOW  process.  A  collection  of  post-processing  utilities  produced  optical 
quality  metrics  to  characterize  the  effects  of  flow  treatment  devices  upon  aero-optical  transmission  quality. 
As  noted  in  Reference  [21],  the  integration  of  the  spectral/parabolic  paraxial  solver  within  the  CFD  model 
was  relatively  efficient,  and  for  most  cases,  reduced  the  processing  speed  of  the  model  by  only  1-3%.  The 
solver  compared  very  well  with  theory /analytical  solutions  where  applicable,  and  demonstrated  spectral 
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convergence  on  certain  problems  [21].  Due  to  the  tight  integration  of  the  paraxial  solver  within  the  CFD 
(OVERFLOW)  model,  some  limitations  and  restrictions  were  encountered.  The  most  restrictive  limitation 
was  the  requirement  that  the  paraxial  equations  be  solved  solely  on  separately  identified  grid  zones  that 
are  part  of  the  overall  overset  grid  system.  This  means  that  the  beam  propagation  directions  and  spans 
have  to  be  determined  during  the  mesh  generation  phase,  and  cannot  be  changed  afterwards.  In  addition, 
during  the  grid  compositing  process,  the  beam  path  grids  cannot  be  interpolated  from  any  other  donor  grids 
placing  unnecessary  restrictions  on  grid  topology  and  compositing.  The  procedures-based,  tightly  integrated 
approach  within  the  CFD  model  resulted  in  a  standalone  model  that  was  difficult  to  embed  in  other  models, 
such  as  adaptive-optics  or  propagation  models. 

In  order  to  provide  a  cleaner  computational  interface  to  other  CFD  models,  we  have  developed  an  Aero- 
Optics  software  class  library,  and  have  used  it  to  embed  the  spectral/parabolic  paraxial  beam  solver  within 
both  an  overset,  structured  solver  (OVERFLOW)  and  an  unstructured,  finite- volume  solver  (AVUS).  The 
integration  within  the  CFD  models  is  made  as  straightforward  as  possible,  by  providing  four  Application 
Programming  Interfaces  (API)  called  by  the  CFD  codes.  A  mesh  pre-processing  phase  finds  interpolation 
stencils  and  their  weights  for  a  given  mesh  system  and  beam  grid  layout,  and  stores  this  data  in  a  commonly 
accessed  file.  To  speed  up  the  software  development,  this  common  file  format  has  been  constructed  using 
a  specialization  of  the  NetCDF  library,  and  provides  direct  file  access  of  commonly  shared  data  by  all  of 
the  applications  in  the  aero-optics  library  framework.  To  provide  maximum  flexibility,  the  paraxial  beam 
solver  may  be  operated  in  a  standalone  fashion,  where  it  uses  the  CFD  generated  density  field,  stored  only 
at  the  limited  number  of  points  needed  to  interpolate  the  data  from  the  CFD  mesh  to  the  beam  grids. 
Post  processing  tools  are  provided  that  perform  wave  front  processing  tasks,  and  are  used  to  characterize  the 
aero-optical  fields  in  terms  of  metrics,  such  as  beam  jitter  and  bore  sight  error.  Other  tools  are  provided  that 
compute  wave  front  error  using  Zernike  polynomials,  and  may  also  be  used  to  characterize  the  aero-optical 
distortions. 

In  the  sections  that  follow,  we  describe  the  spectral/parabolic  operator  splitting  method  used  to  solve 
the  parabolized  Maxwells  equations.  The  top-level  framework  of  the  library  is  described  next,  showing  the 
typical  data  paths  and  steps  taken  to  produce  beams  and  post-process  them.  A  detailed  description  of  the 
mesh  pre-processor  is  then  given,  followed  by  a  description  of  the  CFD  model  API  provided  by  the  aero- 
optics  library  and  how  we  used  these  to  integrate  OVERFLOW  and  AVUS  with  the  library.  Post  processing 
tools,  including  the  standalone  beam  solver  and  aero-optics  metrics  generation  are  shown  next.  Finally,  we 
demonstrate  the  model  for  the  DES  simulations  past  a  circular  cylinder  using  both  OVERFLOW  and  AVUS, 
and  compare  measured  to  computed  beam  metrics  and  aerodynamic  quantities  for  the  flow  past  a  1-foot 
diameter,  conformal,  hemispherical  turret  mounted  in  the  Air  Force  Academy  wind  tunnel. 


3  Spectral/Parabolic  Solution  of  the  Paraxial  Beam  Equation 


The  aero-optics  model  is  based  upon  solving  the  parabolized  Helmholtz  equation,  also  referred  to  as  the 
paraxial  beam  equation  [17,  22].  This  is  a  reduced  form  of  Maxwells  equations  which  govern  electromagnetic 
propagation  through  general  media.  For  the  electromagnetic  propagation  considered  here,  the  medium,  air, 
is  considered  to  be  non-conducting  with  a  constant  magnetic  permeability  [3].  In  addition,  since  the  time  scale 
for  the  electromagnetic  propagation  is  orders  of  magnitude  faster  than  the  flow  time  scales,  the  temporal 
variation  of  the  electromagnetic  field  is  assumed  to  be  negligible.  Depolarization  is  also  neglected  [17,  22], 
since  the  propagating  wavelength  is  assumed  to  be  much  smaller  than  the  turbulence  length  scale  [17]. 
Neglecting  propagation  normal  to  the  beam  direction,  assuming  planar  wave  propagation  and  applying  the 
paraxial  approximation  results  in  the  following,  where  the  coordinates  are  scaled  by  the  propagation  beam 
wavelength  (£  =  x / A,  77  =  y / A,  £  =  z/ A): 
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The  Gladstone-Dale  constant  relates  the  index  of  refraction  to  the  fluid  density  for  a  given  wavelength  (in 
meters)  as 

n  =  1  +  Kgdp  (2) 
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where 


Kgd  =  2.24  x  10“4 
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We  solve  the  paraxial  beam  equation  (1)  by  a  spectral/parabolic  method  using  operator  splitting [21].  The 
first  parabolic  operator  is  termed  a  phase  shift  advance  since  it  only  alters  the  phase  of  the  planar  wave, 
which  is  followed  by  the  Fourier-based  wave  advance,  which  represents  the  Laplacian  operator  in  the  beam 
tangential  direction  via  a  Fourier  decomposition,  and  advances  the  spectral  coefficients  in  the  propagation 
direction.  The  spectral/parabolic  method  marches  the  optical  beam  through  the  density  field  by  solving 
successively  a  phase  shift  (4)  and  a  wave  advance  (5)  operator. 
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The  final  form  of  the  wave  advance  equation  in  spectral  space  is: 

Akl(C  +  AC)  = 

The  phase  shift  equation  is  advanced  in  physical  space  as: 

Akl((  +  A()  =  Akl(0ej”[(^2-1lAC 
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The  paraxial  solver  has  been  validated  previously  by  comparison  to  accepted  propagation  code  results  for  a 
turbulent  field,  as  well  as  to  theory  for  propagation  of  a  beam  through  a  vacuum.  A  Gaussian  beam  may  be 
expressed  in  non-dimensional  form  as: 
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and  r~o  =  ro/A  and  /  =  //A.  We  compute  the  solution  on  spectral  grids  with  densities  of  (Nx,Ny)  ranging 
from  8  to  1024.  The  Gaussian  beam  parameters  are  Ao  =  1,  /  =  20,  000,  r~o  =  10,  A  =  0.1m,  Lx  =  Ly  =  400 
and  Lz  =  1000.  Error  norms  are  computed  that  compare  the  magnitudes  of  the  computed  and  exact  solution 
at  the  propagation  distance  noted  above,  and  are  shown  in  Figure  1  plotted  against  the  mesh  size  h  =  Lx/Nx. 
As  can  be  seen  from  the  plot,  the  error  norms  exhibit  spectral  convergence.  As  defined  in  Boyd  [23],  spectral 
convergence  (which  is  synonymous  to  exponential  convergence)  results  in  an  error  decreasing  faster  than 
e  ~  1/7V^  for  any  power  k. 

We  note  that  although  we  have  taken  an  operator  splitting  approach  here  to  solve  the  paraxial  equation, 
it  is  a  simple  matter  to  replace  the  operator  splitting  approach  with  an  unsplit  method,  and  use  the  spectral 
model  to  represent  the  beam  variation  in  the  plane  of  propagation,  and  use,  say,  a  Runge-Kutta  or  other 
explicit  method  to  march  in  the  propagation  direction.  An  adaptive  Runge-Kutta  method  that  uses  a 
compact  scheme  has  been  successfully  applied  for  solving  the  paraxial  beam  equation  [24]. 


4  Aero-Optics  Class  Library,  Framework  and  Data  Flow 

The  Aero-Optics  capability  is  built  upon  a  series  of  public-domain  libraries  and  specially  written  classes 
which  form  a  software  library.  Using  this  library,  a  number  of  standalone  applications  have  been  written, 
including  the  API  used  by  the  CFD  models  to  invoke  the  software  library.  The  block  diagram  shown  in 
Figure  3.1  outlines  the  different  components. 
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Gaussian  Beam,  Wave  Advance  Only,  Lx=Ly=40 

Grid  Refinement  Study  with  paraxialExecutorClass 


Figure  1:  Convergence  rate  for  the  Gaussian  beam  problem. 
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Post  Processing 
Tools 

preprocessor 

cfdPublicAPI 

libAO:  C++  Class  Library  for  Aero-Optics 

dfi file Jo:  NetCDF  based  common  file  format  and  library 

Figure  2:  Aero-Optics  Class  Library  Components  and  Framework 

4.1  Library  Components 

Referring  to  Fig.  2,  the  components  of  the  Aero-Optics  capability  are  (from  bottom  to  top): 

•  df_file_io:  This  is  a  C++  library  and  namespace  which  uses  a  specialization  of  the  NetCDF  [25]  common 
file  format  and  library  to  provide  a  common  set  of  API  that  are  used  to  read/write  named  data 
from/to  a  commonly  accessed  file.  The  data  supported  include  vectors  (singly  subscripted  arrays)  of 
basic  data  types  (unsigned  ints,  ints,  floats,  doubles)  as  well  as  complex  data  and  strings.  Access  by 
Fortran  applications  is  also  provided.  By  allowing  direct  access  to  data  in  a  commonly  shared  file,  the 
software  development  process  has  been  sped  up  considerably.  Relevant  classes  in  the  libAO  library 
have  readFromCommonFile  and  writeToCommonFile  functions. 

•  libAO:  This  is  the  C++  class  library  that  contains  all  the  classes  and  functions  used  by  the  Aero-Optics 
applications.  The  standalone  applications  (described  below)  are  all  constructed  from  classes  supplied 
by  this  library.  The  paraxial  solver  invoked  by  the  CFD  models,  and  all  the  ancillary  processing,  is 
also  handled  by  classes  that  make  up  this  library.  A  sampling  of  some  of  the  classes  includes: 

—  paraxialExecutorClass:  This  class  controls  and  marshals  the  individual  paraxial  solvers,  and  is 
invoked  by  the  standalone  paraxial  solver  and  the  CFD  model  API. 

—  paraxialSolverClass:  This  class  encapsulates  the  paraxial  solver  itself. 

—  boxGridClass::beamGridClass:paraxInterpolatorClass:  This  is  a  set  of  derived  classes  which  are 
used  to  define  beam  grid  extents  and  provide  a  fast  interpolator  to  be  used  by  the  parabolic  solver 
to  find  the  density  at  the  spectral  grid  points  as  the  wave  is  propagated  through  the  beam  grid. 
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•  preprocessor:  This  performs  mesh  pre-processing  for  given  beam  propagation  paths  and  interpolator 
grid  densities  to  find  donor  interpolants  from  the  CFD  mesh,  and  store  them  in  the  commonly  accessed 
file.  Section  5  describes  the  pre-processing  steps  and  approach. 

•  cfdPublicAPI:  Four  simple  functions  are  provided  that  are  invoked  by  the  CFD  model.  The  functions 
are  described  fully  in  Section  6,  and  permit  invocation  of  the  paraxial  beam  propagator,  and  store 
donor  density  data  at  the  limited  set  of  unique  CFD  donor  points  used  to  define  the  density  on  the 
beam  propagator  grids. 

•  Standalone  Paraxial  Solver:  The  paraxial  beam  propagator  can  be  run  in  a  standalone  manner,  and 
uses  the  donor  density  data  stored  when  the  CFD  model  is  run.  The  input  to  the  standalone  paraxial 
solver  is  identical  to  that  when  it  is  invoked  from  within  the  CFD  model,  and  produces  identical  data 
files  as  well. 

•  Post  Processing  Tools:  A  variety  of  post  processing  tools  are  provided  that  process  propagated  beams 
and  density  fields  and  produce  aero-optical  metrics,  waveforms  and  related  information.  The  data 
generated  by  these  tools  includes: 

—  Wavefront  Errors:  Zernike  polynomials  are  used  to  remove  piston,  tip  and  tilt  from  the  beams,  and 
compute  unsteady  wave  front  error  wave  forms,  that  are  further  processed  to  compute  statistical 
descriptions. 

—  Aero-Optical  Metrics:  Quantities  such  as  jitter,  bore  sight  error  and  OPD  are  computed  by 
processing  the  beams  in  time. 

In  addition  to  the  libraries  noted  above,  we  make  use  of  a  number  of  public-domain  libraries  as  well, 
including: 

•  fftw3:  The  FFTW3  library  [26]  is  used  to  perform  the  Discrete  Fourier  Transformation  (DFT)  between 
physical  and  spectral  space.  It  is  publicly  available  and  provides  a  very  useful  set  of  API  for  performing 
DFT  related  operations. 

•  boost:  The  boost  C++  library  (http://www.boost.org/)  is  a  widely  used  set  of  classes  and  libraries 
and  contains  a  variety  of  very  useful  classes  that  are  an  addition  to  the  standard  namespace  C++ 
libraries. 

•  stdmamespace:  The  std  namespace  distributed  with  the  standard  C++  language  maintains  some  very 
useful  container  classes  and  related  functionality. 

4.2  Data  Processing  Work  Flow 

The  data  flow  for  the  Aero-Optics  capability  is  based  upon  a  number  of  shared,  commonly  accessed  files  that 
are  produced  at  various  phases  of  the  processing.  The  typical  processing  steps  are  as  follows: 

1.  Mesh  Pre-Processing:  Given  a  CFD  grid  and  definition  of  the  beam  interpolator  grids  (described  in 
Section  5),  the  mesh  pre-processor  (preprocessor)  finds  the  unique  set  of  CFD  donor  data  points  and 
all  the  interpolants  and  weights  to  interpolate  the  CFD  solution  onto  the  set  of  beam  interpolator 
grids.  This  data  is  stored,  for  the  N-th  beam  grid,  in  a  data  file  (beamGrid.interp.N.nc) 

2.  CFD  Phase:  In  this  phase,  the  CFD  model  is  run.  The  CFD  code  input  is  modified  by  the  addition 
of  a  small  namelist  (for  the  OVERFLOW  and  AVUS  codes),  which  tells  the  solver  when  to  invoke  the 
beam  propagator  (if  selected),  and  at  what  frequency  to  call  it,  as  well  as  the  name  of  the  common  file 
which  defines  the  propagator  input  data.  The  propagator  input  data  file  defines  the  number  of  paraxial 
solvers  to  run,  and  associates  the  (pre  processed)  beam  grid  interpolator  files  with  the  paraxial  solvers. 
This  approach  allows  the  re-use  of  beam  grids  for  different  solvers,  so  that,  for  example,  different 
waveforms  can  be  run  through  the  same  propagation  path.  The  solver  writes  out,  for  each  beam, 
the  CFD  donor  density  values,  for  each  time  step,  to  a  data  file  (rhoDonor.N.nc).  If  the  propagator 
is  turned  on,  the  propagated  beam,  for  each  beam  grid  ,  for  each  time  step,  is  written  to  a  data  file 
(paxaxSolver.N.beam).  These  two  CFD  model  produced  files  (rhoDonor.N.nc  and  paraxSolver.N.beam) 
are  used  subsequently  by  the  following  standalone  processors  and  post-processors. 
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3.  (Optional)  Standalone  Paraxial  Solver:  A  standalone  paraxial  solver  is  available  that  reads  a  propa¬ 
gator  input  file  (that  defines  the  beam  grids  and  data  needed  by  each  paraxial  solver),  and  for  each 
beam  grid,  reads  and  processes  the  donor  density  data,  propagates  the  given  beam  waveform  for  the 
prescribed  time  steps,  and  deposits  the  propagated  beams  into  the  beam  file  (paraxSolver.N.beam). 
This  allows  one  to  re-use  the  computationally  expensive  CFD  data. 

4.  Post  Processing  Phase:  The  post  processors  (described  in  Section  8)  are  used  to  process  given  beam 
files  (ie;  propagated  waveforms)  and  produce  aero-optical  quality  metrics  and  descriptions.  The  post 
processors  operate  on  given  propagator  input  files  (again,  defining  beam  grids,  interpolator  data  and 
rhoDonor  data),  and  produce  the  metrics. 

The  following  sections  describe  the  different  components  in  more  detail. 


5  Mesh  Pre-Processor 

In  order  to  make  the  paraxial  solver  computationally  efficient,  it  must  be  able  to  quickly  compute  the  index 
of  refraction  on  the  integration  points  of  the  spectral  grid  planes  as  they  are  marched  in  the  propagation 
direction.  The  index  of  refraction  is  directly  related  to  the  density  via  the  Gladstone-Dale  constant,  so  this 
means  that  a  fast  interpolation  method  is  needed  to  interpolate  the  density  from  the  CFD  grid  onto  the 
spectral  grid  as  it  is  marched  in  the  propagation  direction.  This  is  accomplished  by  a  pre-processing  step  that 
finds  the  CFD  mesh  donor  interpolants  (donor  cells  and  weights)  on  a  Cartesian  grid  that  lies  in  the  beam 
propagation  path.  This  Cartesian  grid  is  aligned  in  space  with  what  is  termed  a  beam  grid.  The  Cartesian 
grid  itself,  upon  which  the  CFD  donor  data  and  interpolants  are  precomputed  is  termed  an  interpolator 
grid.  The  paraxial  solver  marches  in  the  beam  propagation  direction,  and  needs  the  index  of  refraction  at 
the  points  that  are  termed  the  spectral  grids.  The  following  diagrams  explain  the  relationship  between  beam 
grids,  interpolator  grids  and  spectral  grids. 


Figure  3:  Illustration  of  beam  grids,  interpolator  grids  and  spectral  grids 

The  preprocessor  processes  a  given  set  of  interpolator  grids  for  a  particular  CFD  mesh  and  finds,  for  every 
point  on  the  interpolator  grid,  the  nearest  CFD  grid  element  (nodes  for  structured  grids,  cells  for  unstructured 
grids) .  For  each  of  these  elements  a  corresponding  stencil  is  found  dependent  upon  the  particular  flow  solver 
type  (ie;  structured  overset,  or  unstructured),  and  the  weights  needed  to  interpolate  the  CFD  solution  onto 
the  interpolator  grid  point  are  computed.  Once  all  the  donor  points  (ie:  CFD  grid  points  that  participate 
in  the  interpolation)  are  found,  a  set  of  unique  donor  points  is  created.  This  set  of  unique  donor  points  is 
what  needs  to  be  gathered  by  the  particular  CFD  solver,  and  is  passed  to  the  paraxial  beam  model  during 
the  time-varying  CFD  computations.  Since  this  set  of  unique  donor  points  is  a  significantly  smaller  subset 
of  the  entire  CFD  grid,  we  provide  the  ability  to  store  the  density  at  these  points  during  the  transient  CFD 
solution,  which  form  the  rhoDonor. N.nc  files.  The  pre-processor  uses  a  modification  of  an  ADT  tree [27]  that 
was  supplied  to  us  by  Ralph  Noack  at  PSU.  This  space  partitioning  tree  permits  very  fast  nearest  neighbor 
searches  that  are  needed  to  find  the  donor  grid  points  for  the  interpolator. 
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For  either  structured  overset  or  unstructured  grids,  the  interpolation  from  the  CFD  grid  to  the  interpo¬ 
lator  grid  is  represented  as  a  sum: 


N 


Pi  ^  ^  &  Pn 

n=l 


(10) 


Structured  overset  grid  weights  are  found  using  a  hexahedral,  iso-parametric  representation  in  local,  natural 
coordinates.  For  a  given  interpolator  point,  the  natural  coordinates  are  found  by  solving,  via  Newton 
iteration: 


X  =  (£,??,  C )xn 

n 

y  =  '52(t>n  (£,??,  C )Vn 

n 

Z  =  'Yh(t)n  (t,y,()Zn 

n 

Then,  the  interpolation  weights  are  readily  found  from  the  iso-parametric  interpolation  basis  functions. 

Unstructured  grid  interpolator  weights  are  found  by  first  finding  the  grid  nodes  that  make  up  the  given 
nearest  cell  to  the  interpolator  point  (which  is  found  by  the  ADT  nearest  neighbor  search)  using  cell-to-node 
connectivity.  This  cloud  of  points  is  then  used  in  a  linearity  preserving  reconstruction  scheme  to  find  the  grid 
node  weights  in  the  donor  value  reconstruction.  For  the  interpolations  used  here,  we  find  the  weights  using 
the  linearity  preserving  Laplacian  scheme  [28]  which  finds  the  weights  so  that  an  interpolation  operator, 

£(/)  =  (/-/„)  =  0  (12) 

n 

is  exactly  satisfied  when  /  is  a  linear  function  of  the  coordinates.  That  is,  find  the  Lagrange  multipliers  in  is 
exactly  satisfied  when  f  is  a  linear  function  of  the  coordinates.  That  is,  find  the  Lagrange  multipliers  An  in 

U}n  =  1  +  Xx  (x  -  xn)  +  A y(y-  yn)  +  \z{z-  zn)  (13) 

so  that 

L(x)  =  L(y )  =  L(z)  =  0  (14) 

The  resulting  3x3  system  of  equations  is  inverted,  which  then  gives  the  Lagrange  multipliers,  which  then 
yield  the  interpolator  weights. 

Figure  4  shows  an  example  CFD  grid  with  an  interpolator  grid,  and  the  interpolated  solution  overlayed 
with  the  CFD  grid  solution  for  an  OVERFLOW  and  AVUS  solution  for  the  unsteady  flow  over  a  circular 
cylinder. 


(ii) 


6  Application  Programming  Interface  to  CFD  Models 

The  API  supplied  by  the  Aero- Optics  library  to  be  called  by  the  CFD  models  has  been  made  as  minimal  as 
possible.  There  are  a  total  of  four  API  called  by  CFD  model,  and  are  represented  by  the  following  pseudo¬ 
code: 

initializer  araxialB  earn  Solvers  ( commonFormatFileN  ame  ); 
for  each  timestep  in  timesteps 
for  each  beamNum  in  beams 

setRho  Donor  (beamNum, rho  Donor,  timestep) ; 
if(paraxSolve )  runParaxialB  earn  Solver  (beamNum,  timestep ); 
finalizeParaxialBeamSolvers  (); 


•  Initializes eamSolvers  instantiates  the  paraxial  solver  executor  and  instructs  it  to  load  all  its  data 
from  the  supplied  common  file  named  commonFormatFileName. 
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(a)  (b)  (c) 


Figure  4:  Streamwise  aligned  cut  through  the  CFD  grid  and  interpolator  grid  (left),  OVERFLOW  solution 
and  interpolated  density  (middle)  and  AVUS  solution  and  interpolated  density  (right)  at  a  particular  time 
step  of  an  unsteady  DES  calculation 


•  setRhoDonor  transfers  the  CFD  supplied  donor  density  values  to  the  interpolator  grid,  which  then 
loads  up  the  interpolator  with  the  density  data  using  the  precomputed  interpolants.  The  donor  density 
data  is  appended  to  a  local  NetCDF  file 

•  runParaxialB  earn  Solver  runs  the  paraxial  beam  solver  at  this  time  step  for  this  beam  using  pa¬ 
rameters  specified  in  the  common  file. 

•  finalizeP  araxialB  earn  Solvers  performs  all  closeout /finalization  for  the  solvers. 

We  have  integrated  the  paraxial  solver  capability  within  both  OVERFLOW  and  AVUS.  Both  solvers  are 
modified  to  read  in  a  small  namelist  that  defines  whether  or  not  the  paraxial  model  is  to  be  invoked,  whether 
or  not  to  run  the  beam  propagator,  how  often  to  invoke  the  model,  and  the  name  of  the  paraxial  model 
common  file.  A  Fortran90  module  is  created  for  each  solver  to  store  this  data,  as  well  as  the  donor  mesh 
connectivity  data  that  is  read  from  the  pre-processor  produced  files.  Specific  to  each  solver,  a  density  gather 
operation  has  been  written  to  gather  and  store  only  the  donor  density  data,  which  is  passed  to  the  model 
through  the  setRhoDonor  API.  For  the  AVUS  code,  a  specialization  of  the  procedure  that  interpolates  data 
from  cell  centers  to  nodes  has  been  made  to  perform  this  density  gather.  Integration  within  both  CFD 
models  has  been  straightforward,  and  lends  itself  to  integration  in  other  models  as  needed. 

7  Standalone  Paraxial  Beam  Model 

The  paraxial  library  software  framework  has  allowed  the  straightforward  development  of  a  standalone  parax¬ 
ial  beam  solver.  The  standalone  paraxial  beam  solver  is  designed  to  use  data  generated  from  the  CFD  model 
in  the  form  of  the  donor  density  files,  and  produces  identical  beam  files  that  are  the  beams  propagated 
through  the  unsteady  flow  field.  By  following  the  same  data  flow  model,  it  is  straightforward  to  reuse  the 
computationally  expensive,  transient,  CFD  data  for  other  applications  as  well.  The  potential  applications 
include: 

•  Integration  within  Far-Field  Models:  We  see  this  as  being  directly  integrated  within  engineering  level 
models  which  may  be  used  to  evaluate  the  effect  of  the  near  field  aero-optics  upon  the  full  scale 
engineering  system.  Potential  models  include  the  WaveTrain  model  produced  by  MZA  Associates. 

•  Multiple-Path  Propagation:  The  standalone  model  allows  one  to  propagate  a  beam  out  of  the  near 
field,  perform  a  modification  to  this  beam  (through,  for  instance,  an  engineering  model  as  noted  above) 
and  then  propagate  it  back  through  the  near-field,  and  measure  its  beam  quality. 

•  Multiple  Wave  Forms  and  Wavelengths:  The  standalone  model  allows  many  different  wave  forms  at 
different  wavelengths  to  be  modeled  through  the  same  CFD  generated  flow  field. 
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•  Evaluation  of  Beam  Propagator  Methods:  It  is  a  relatively  simple  matter  to  replace  the  factored  parax¬ 
ial  beam  propagator  developed  here  with  an  unsplit  method,  such  as  an  adaptive  RK  method,  similar 
to  that  shown  in  Reference  [24] ,  and  use  the  standalone  approach  to  evaluate  and  compare  the  different 
solution  strategies. 

8  Post-Processing  Tools  and  Models 

A  variety  of  post-processing  tools  have  been  developed  using  the  software  library,  and  are  designed  to  post¬ 
process  given  waveforms  and  density  fields  to  produce  aero-optical  quality  metrics.  For  given  waveforms, 
the  beam  jitter  and  boresight  error  are  constructed  for  each  timestep,  and  output  to  ASCII  files  for  plotting 
and  processing.  In  order  to  compare  with  experimental  data,  wavefront  error  is  computed  using  Zernike 
polynomials  and  statistical  representations  of  the  propagated  beams  are  also  made  available.  As  an  alter¬ 
native  to  using  propagated  beams  to  characterize  the  aero-optical  field,  the  Optical  Path  Difference  (OPD) 
computed  using  a  purely  geometric  integral  approach  is  also  made  available,  and  these  phase  errors  may  also 
be  statistically  processed  by  computing  wavefront  error.  The  sections  below  describe  the  approaches  taken 
for  each  of  the  methods. 

8.1  Jitter  and  Boresight  Error 

Using  beams  contained  in  the  “paraxialSolver.N.beam  files,  beam  jitter  and  boresight  error  are  calculated 
using  beam  power- weighted  integrals  as  follows.  For  each  stored  time  step,  the  beam  centroid  is  found  by 
using  a  power- weighted  area  integral.  For  a  general  function  defined  in  the  imaging  plane,  we  define  a  generic 
power- weighted  area  integral  to  be: 


1(f)  =  [  [  AA*f(x,y)dxdy 

J  x  J  y 


(15) 


where  A  A*  is  the  beam  power  found  by  multiplying  the  beam  magnitude  by  its  complex  conjugate.  Then 
we  find,  at  the  nth  time  slice,  the  centroids  to  be: 


I(x)  I(y ) 

Xn  =  7(1)  ’  Vn  =  7(1) 

The  average  beam  centroid  is  found  from  the  beam  centroids  and  the  N  time  intervals  as: 

N  ,  N 


-  n  FA’  Vn~  jv  F 


Vn 


(16) 


(17) 


n= 1  n=  1 

We  define  jitter  as  the  root  mean  square  of  the  time- varying  beam  centroid  about  the  average  beam  centroid: 


(7  rp  - 


\\  jy  F  (Xn  Xn )> 


\ 


1  N 

x  F  (yn  - 


n=  1 


The  angular  deviation  over  the  propagation  length  is: 


0T  = ; 


*(t)- 


(18) 


(19) 


We  define  boresight  error  using  the  difference  between  the  beam  centroid  and  image  center.  The  translational 
error  is: 

^ x,n  =  %n  ^0, image-)  ^y,n  =  Vn  VO,  image  (20) 

where  the  boresight  error  is 

\J  ^x,r 


+  e2 

:,n  1  y,n 


(21) 
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8.2  Wavefront  Error  Processing 

We  provide  the  ability  to  compute  the  wavefront  error  for  each  beam,  and  then  to  statistically  process 
these  wavefront  errors.  The  approach  taken  here  to  produce  the  wavefront  error  is  to  compute  the  phase 
of  the  beam,  and  then  using  Zernike  polynomials  compute  and  subtract  from  the  phase  the  piston,  tip  and 
tilt.  Zernike  polynomials  are  a  set  of  orthogonal  polynomials  defined  on  a  unit  disk  that  are  often  used  in 
optics  to  characterize  waveforms  in  a  consistent  manner,  of  which  the  piston,  tip  and  tilt  are  the  first  three 
polynomials.  We  use  ANSI  standard  Zernike  defined  for  the  n-th  radial  and  m-th  azimuthal  order  as: 


(p)  cos  (mO) ;  m  >  0  1 

— (p)  sin  (mO) ;  m  <  0  J 


(22) 


where 


and 


(n—\m\)/2 

=  E 


(-!)>-«)! 


-n—2s 


s= 0 


s\  [(n  +  \m\)  /2  —  5!]  \{n  —  \m\)  /2  —  s\] 


NN  = 


2n  T  2 


1  +  ^mO 

For  a  given  wave  form,  the  wavefront  error  is  computed  by  the  following: 


(23) 


(24) 


1.  Compute  a  no  flow  beam  phase  (0^)  by  propagating  the  beam  through  the  specified  propagation 
distance  using  a  constant  density  field,  set  at  the  freestream  density  value. 

2.  Compute  the  difference  between  the  n-th  timesteps  beam  phase  and  the  no  flow  phase:  0i  =  0  —  q 

3.  Compute  the  piston  (0P),  tip  {(j>tiP)  and  tilt  (4>tut)  using  the  orthogonality  properties  of  the  Zernike 
polynomials,  and  subtract  these  from  the  phase  difference:  02  =  0i  —  (0P  +  4>tip  +  4>tilt) 

4.  Find  the  mean  phase  error:  02  =  (02) 

5.  Find  the  n-th  timesteps  wave  front  error  by  subtracting  the  mean  from  the  phase  error:  0  =  02  —  02 


This  wave  front  error  is  then  statistically  processed  by  computing  a  mean  and  the  root  mean  square  deviation 
from  the  mean.  The  procedure  described  above  requires  the  definition  of  an  aperture  mask,  since  the  Zernike 
polynomials  are  defined  on  a  unit  radius  disk. 


9  Optical  Path  Difference-based  Wavefront  Error  Processing 

Here,  instead  of  propagating  a  beam  through  the  beam  grid,  we  integrate  the  density  along  the  propagation 
direction,  generating  the  Optical  Path  Difference  (OPD).  The  OPD  is  defined  as: 

[Lz 

OPD  =  Kgd A  /  pd(  (25) 

Jo 

The  OPD  is  related  to  the  phase  shift  since 

r^z  97 r 

A 0  «  2t rKgd  /  pd(  =  -—OPD  (26) 

Jo  * 

The  same  wavefront  error  processing  shown  above  is  used  to  compute  a  statistical  description  of  the  OPD- 
based  wavefront  error. 
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Figure  5:  Unsteady  flow  past  a  circular  cylinder.  Entropy  measure  contours  (left),  Mach  contours  in  a 
spanwise  plane  (middle),  locations  and  directions  of  beam  propagator  grids  (right). 


10  Demonstration:  Unsteady  Flow  Past  a  Circular  Cylinder 

The  paraxial  solver  is  first  demonstrated  for  the  unsteady  flow  past  a  circular  cylinder  using  both  the 
OVERFLOW  and  AVUS  CFD  models.  The  cylinder  diameter  is  2  inches  and  the  free  stream  Mach  number 
is  0.2  with  static  conditions  corresponding  to  1  atmosphere  at  288  degrees  K.  The  time  step  is  chosen  to  be 
At  =  8.5  x  10  6  seconds  which  resolves  one  shedding  period  with  200  time  steps  when  using  an  estimate  of 
the  shedding  frequency  based  upon  a  Strouhal  number  of  0.22.  The  cylinder  span  is  taken  to  be  10  inches. 
Both  models  use  the  same  mesh,  which  is  a  structured,  single  domain  grid  of  101  x  151  x  201  grid  points 
in  the  spanwise,  radial  and  azimuthal  directions.  For  the  AVUS  code,  a  mesh  processor  was  written  that 
converts  the  structured  grid  into  AVUS  input  format.  The  mesh  normal  to  the  wall  is  clustered  to  give  the 
first  cell  center  a  value  in  log-normal  coordinates  of  y+  ~  30  which  is  adequate  for  both  models  using  wall 
functions.  Both  codes  were  run  with  the  Spalart-Allmaras  Detached  Eddy  Simulation  [29]  turbulence  model 
for  16,000  time  steps,  followed  by  1,000  time  steps  with  the  paraxial  model  activated.  OVERFLOW  was 
run  using  a  dual  time  stepping  scheme,  the  HLLC  approximate  Riemann  solver  and  third-order  upwinding 
with  the  default  limiter.  AVUS  was  run  with  dual  time  stepping,  a  least-squares  reconstruction  scheme,  the 
HLLC  approximate  Riemann  solver  and  limiting  turned  off.  The  paraxial  model  recorded  the  donor  density 
at  5  time  step  intervals,  which  was  subsequently  used  by  the  post  processing  models  to  compute  various 
aero-optical  beams  and  beam  metrics.  The  beam  grids  shown  on  the  right  of  Fig.  5  shows  that  the  first  two 
grids  are  oriented  so  that  the  beam  shines  upwards,  through  the  shear  layer  sheet,  while  the  third  beam 
grid  shines  along  the  spanwise  direction,  just  aft  of  the  cylinder.  The  interpolator  grids  associated  with  each 
beam  grid  were  61  x  61  x  61,  from  which  the  pre-processing  found  63,895,  15,440  and  277,419  donor  points. 
The  pre-processor  took  71  and  99  seconds  to  process  the  overset  and  unstructured  grid  systems,  respectively, 
to  find  the  donor  points  and  interpolator  weights. 

Using  the  donor  density  deposited  every  5  time  steps,  top  hat  beam  profiles  of  roughly  5  inch  diameter 
were  propagated  through  the  beam  grids  for  wavelength  A  =  1  /im  on  spectral  grids  of  256  x  256  using  the 
standalone  propagator.  Figure  6  shows  the  wave  front  errors  computed  from  the  propagated  beams  using 
the  procedure  outlined  in  Section  8.2  for  a  mask  with  a  diameter  of  4.5  inches. 

An  example  of  post  processing  the  beams  for  aero-optical  propagation  metrics  is  shown  in  Fig.  7,  where 
we  compare  the  computed  jitter  for  the  three  beam  grids  for  the  OVERFLOW  and  AVUS  computations.  It 
is  interesting  to  note  that  the  OVERFLOW  model  gives  results  that  appear  to  have  a  finer  turbulent  scale 
than  AVUS.  This  is  evident  from  the  jitter  shown  in  Fig.  7,  where  the  levels  of  jitter  from  OVERFLOW 
are  higher  than  those  from  AVUS.  Examining  the  wave  front  error  from  beam  grid  1  in  Fig.  6  shows  that 
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(d)  (e)  (f) 


Figure  6:  Comparison  of  Wave  Front  Errors  for  OVERFLOW  (top)  and  AVUS  (bottom)  on  beam  grids  1, 
2  and  3  (across).  Shown  are  computed  rms  values  of  the  wavefront  error  in  microns 


Figure  7:  Comparison  of  jitter  for  the  three  beam  grids  between  OVERFLOW  and  AVUS. 
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(a)  (b) 


Figure  8:  Comparison  of  Mach  contours  mid-span  for  OVERFLOW  (left)  and  AVUS  (right).  Both  contours 
made  with  same  range  and  levels  of  Mach  number  (21  levels  from  0  to  0.2). 


both  codes  predict  spanwise  turbulent  structures,  but  those  from  the  AVUS  produce  higher  wave  front  error 
near  the  cylinder  than  the  OVERFLOW  model.  These  get  dissipated  more  quickly  downstream  at  beam 
grid  2  from  AVUS  than  from  OVERFLOW,  while  the  third  beam  grid,  which  shines  along  a  spanwise  beam 
next  to  the  cylinder,  shows  a  much  more  compact /necked  structure  from  AVUS.  We  speculate  that  this  may 
be  caused  by  differences  in  the  implementation  of  the  HLLC  scheme,  as  well  as  differences  in  the  numerics 
related  to  the  turbulence  models.  Shown  in  Fig.  8  are  Mach  contours  mid-span  of  the  cylinder  at  the  same 
time  step  from  OVERFLOW  and  AVUS.  The  AVUS  results  show  a  more  regular  vortex  street  than  the 
OVERFLOW  results. 


11  Validation:  1  Foot  Conformal,  Hemispherical  Turret 

The  combined  aero-optics  CFD  modeling  capability  has  been  validated  by  comparison  to  both  aerodynamic 
and  aero-optic  data  taken  in  a  wind  tunnel  test  in  the  Air  Force  Academy  3ft.  x  3ft.  subsonic  wind 
tunnel  [13,  30].  A  1  foot  diameter  hemispherical  turret  was  mounted  in  this  tunnel  on  a  4.5  inch  high 
cylinder,  which  itself  is  mounted  to  a  wall  of  the  wind  tunnel.  Details  behind  the  experiment  can  be  found  in 
References  [13,  19,  30],  which  has  proven  the  basis  for  a  number  of  CFD  model  validation  efforts,  including 
those  shown  in  References  [31,  14]  and  others.  The  data  provided  from  the  experiment  includes: 

•  Aerodynamic  Data: 

—  Hot  Wire  Measurements:  Time-averaged  streamwise  velocity  and  streamwise  velocity  fluctuations 
at  five  traverse  locations:  One  upstream  of  the  turret  and  4  downstream,  in  the  wake. 

—  Surface  Pressure  Measurements:  Pressure  coefficients  along  the  hemisphere  surface  on  a  curve 
aligned  with  the  free  stream. 

•  Aero-Optic: 

—  Malley  Probe 

—  2-D  Shack-Hartmann  Wavefront  sensor 

—  8x8  High-Bandwidth  Wavefront  sensor 

For  the  comparisons  shown  here,  we  focus  upon  the  aerodynamic  data  and  the  data  taken  from  the  2-D 
Shack-Hartmann  Wavefront  sensor  for  the  M  =  0.4  test. 
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(a) 


(b) 

Figure  9:  Grid  system  showing  the  upstream  fetch,  wall  clustering  and  nested  wake  grids. 


(a) 


(b) 


Figure  10:  Beam  grids  and  hot  wire  probe  traverse  locations  in  the  CFD  model 
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An  overset  mesh  system  was  developed  with  six  grids  that  models  the  wind  tunnel  test  section,  including 
the  upstream  fetch  of  the  incoming  boundary  layer  to  the  test  section.  The  mesh  along  the  tunnel  walls 
and  turret  walls  was  clustered  to  have  an  initial  spacing  so  that  the  first  cell  is  located  at  a  y+  value  of 
approximately  1,  with  the  first  four  layers  at  this  same  spacing,  followed  by  a  growth  out  to  y  —  0.2  feet 
with  a  stretching  ratio  of  1.15.  Two,  nested,  Cartesian  grids  are  clustered  over  the  turret,  and  are  shown  in 
Fig.  9.  The  total  mesh  has  approximately  14  x  106  nodes. 

To  compare  to  the  experimental  aero-optic  data,  six  beam  grids  are  constructed  to  he  at  elevation  angles 
of  60,  76,  90,  103,  120  and  132  degrees,  and  are  shown  in  Fig.  10,  along  with  the  locations  of  the  hot 
wire  traverses  and  surface  pressure  measurement  locations.  Each  beam  grid  was  used  to  define  interpolator 
grids  with  dimensions  61  x  61  x  101,  which  were  preprocessed  to  produce  the  donor  connectivity  data  and 
interpolators.  Using  a  modified  Strouhal  number,  the  shedding  frequency  of  the  turret  was  estimated,  and 
based  upon  a  single  period,  the  time  step  was  chosen  so  that  200  steps  are  taken  for  each  period.  OVERFLOW 
was  run  with  the  SST-DES  model,  using  dual-time  stepping,  with  the  third-order  upwinding  option  with 
default  limiter,  without  wall  functions,  for  a  total  of  12,500  steps.  The  collection  of  time-averaged  statistics 
began  at  7,500  steps,  and  the  paraxial  model  was  invoked  during  the  last  1,000  steps,  depositing  donor 
density  every  five  time  steps.  The  time-averaged  data  collected  by  OVERFLOW  during  the  DES  calculation 
is  used  to  compute  the  averaged  data  shown  in  Fig.  11  where  only  the  resolved  turbulence  is  used  to  compare 
to  the  velocity  fluctuations.  The  traverse  locations  are  numbered  so  that  location  1  is  upstream,  3  and  5  are 
located  directly  downstream  of  the  turret,  and  2  and  4  are  offset. 

Inspection  of  the  data  in  Fig.  11  shows  that  the  surface  pressure  on  the  turret  is  adequately  predicted, 
although  the  incoming  turbulent  boundary  layer  is  too  thin,  and  the  turbulence  levels  are  underpredicted. 
As  shown  in  the  data  for  probes  2  and  4,  the  turbulence  levels  off- centerline  are  overpredicted,  with  accom¬ 
panying  lower  momentum  regions  near  the  wall  when  compared  to  the  experimental  data.  Along  the  cylinder 
centerline,  the  opposite  behavior  is  shown,  with  the  CFD  model  under-predicting  the  turbulent  fluctuations, 
and  over-predicting  the  velocity  profiles.  This  could  be  due  to  the  low  levels  of  turbulence  being  imposed  on 
the  inflow  to  the  turret  and  warrants  future  investigation. 

The  collected  donor  density  data  was  then  used  to  propagate  beam  waveforms  through  the  beam  grids, 
on  spectral  grids  of  256  x  256  dimensions.  Top  hat  profiles  of  beams  with  wavelengths  of  A  =  1  fim  were 
propagated  for  each  of  the  200  collected  time  steps,  and  wave  front  errors  computed  and  statistically  an¬ 
alyzed.  In  order  to  compare  the  statistical  data,  both  the  experimental  waveforms  and  the  paraxial  beam 
propagator  generated  waveforms  use  the  processing  noted  in  Section  8.2,  from  which  the  mean  rms  and 
standard  deviations  of  the  rms  of  the  waveforms  are  found  (as  noted  in  Reference  [30])  as: 


OPDrms  =  y/{E(<f)  -£(0)2) 


(27) 


The  computed  wave  front  error  is  used,  and  the  spatial  operator,  E,  is  defined  as: 


E  (/  (r,  6,t)) 


Sr  Jo  frdrffl 

Jr  Jo  rdrd0 


(28) 


Here,  r  varies  over  the  aperture  mask,  taken  to  be  4.5  inches  in  diameter.  This  procedure  is  used  on  the 
experimental  data  contained  in  the  Wave  Front  Sensor  data  files  accompanying  the  experimental  data  set, 
as  well  as  the  beam  produced  phase  differences  found  from  the  beam  propagator  model,  and  produces  data 
similar  shown  in  Reference  [30].  Figure  12  compares  the  measured  and  predicted  statistics  by  plotting  mean 
wave  front  error  root  mean  square  values  with  the  standard  deviation  added  and  subtracted  to  the  mean. 

What  is  noted  from  the  comparison  is: 


•  The  mean  wave  front  errors  appear  to  be  well  predicted,  except  for  the  last  beam  at  132  degrees. 

•  The  mean  wave  front  error  deviation  is  underpredicted  at  the  60  and  76  degree  elevations  (upstream), 
but  well  predicted  at  the  other  elevations. 


This  behavior  can  be  attributed  to  the  underpediction  of  the  turbulence  levels  upstream.  Other  contributing 
factors  could  be  sub-iterative  convergence,  choice  of  upwind  schemes,  limiters  and  turbulence  model  related 
issues.  These  could  be  investigated  in  subsequent  studies. 
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Probe  Location  4.c  Probe  Location  5 


(e)  (f) 


Figure  11:  Centerline  surface  pressure  comparison  (upper  left)  and  averaged  streamwise  velocity  and  stream- 
wise  velocitv  fluctuations  at  the  traverse  locations 
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Figure  12:  Comparison  of  measured  to  predicted  wave  front  distortions.  Shown  are  the  mean  values  +/-  the 
standard  deviations 


12  Conclusions  and  Recommendations 

An  object-oriented  software  library  has  been  developed  that  is  used  to  integrate  a  Paraxial  Beam  solver 
within  both  an  overset  (OVERFLOW)  and  an  unstructured  (AVUS)  CFD  model.  The  Paraxial  Beam 
equation  is  a  parabolized  form  of  Maxwells  equations,  and  is  used  to  propagate  optical  waveforms  through 
the  turbulent,  unsteady  flow  fields  produced  by  the  structured  overset  and  unstructured  solvers.  Previous 
work,  where  a  tightly  coupled,  procedures  based  integration  was  made  with  OVERFLOW,  identified  a  number 
of  deficiencies,  which  have  been  addressed  with  this  work.  To  achieve  this  new  capability,  an  object-oriented, 
software  library  approach  has  been  undertaken  and  has  been  used  to  develop  the  following  components;  a 
standalone  mesh  preprocessor,  an  easily  integrated  Paraxial  Beam  solver,  a  standalone  Paraxial  Beam  solver 
and  an  assortment  of  post-processing  tools.  The  new  capability  has  been  integrated  within  the  OVERFLOW 
and  AVUS  CFD  solvers,  and  demonstrated  for  evaluating  the  Aero-Optical  propagation  through  unsteady, 
turbulent  flows  past  a  circular  cylinder  and  validated  by  comparison  to  aerodynamic  and  aero-optical  data 
for  a  conformal,  hemispherical  turret. 

Based  upon  our  work,  we  recommend  the  following  areas  to  be  investigated: 

•  Evaluate  CFD  Model  Accuracy  and  Resolution  Requirements.  For  the  wavelengths  under  considera¬ 
tion,  it  is  critical  to  capture  the  small  scale  turbulent  structures  in  order  to  provide  accurate  near-field 
propagation  characteristics.  A  compact-scheme  based  CFD  code,  such  as  FDL3DI,  could  be  needed 
to  model  this  behavior  properly,  but  until  a  systematic  study  is  conducted,  this  is  purely  speculator^ 
We  recommend  performing  such  a  study  by  first  integrating  the  paraxial  model  into  FDL3DI,  and 
comparing  canonical  flows  with  various  turbulence  scales  produced  with  varying  levels  of  discretization 
accuracy  in  order  to  ascertain  these  effects.  The  hemispherical  turret  cases  (both  flat  windowedl7  and 
conformall6)  can  provide  experimental  data  on  which  to  base  the  comparisons  as  well. 
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•  Evaluate  Paraxial  Model  Accuracy  and  Resolution  Requirements:  The  effect  of  paraxial  model  resolu¬ 
tion  (marching  step  size  and  spectral  grid  dimensions)  must  be  understood  to  characterize  the  effect  of 
the  paraxial  beam  solution  accuracy  upon  the  predicted  near- field  aero-optical  distortions.  Implemen¬ 
tation  of  a  non-factorized  scheme  should  be  investigated,  and  where  possible,  techniques  such  as  the 
Method  of  Manufactured  Solutions,  as  shown  in  reference  27,  should  be  used  to  formally  characterize 
the  accuracy  of  the  beam  propagation  through  non- vacuum  fields. 

•  Integration  with  Engineering  Scale  Propagation  Models:  The  true  effect  of  near-field  optical  distortions 
can  only  be  measured  in  terms  of  their  impact  upon  the  total  system  performance.  This  is  most  readily 
achievable  by  integrating  the  results  from  near  field  CFD  produced  aero-optical  distortions  into  a  large- 
scale  engineering  propagation  model  which  models  the  entire  end-to-end  process.  The  approach  we  have 
taken  here,  by  storing  donor  density  data  in  commonly  accessed  files,  and  providing  for  standalone 
paraxial  beam  propagators,  forms  the  basis  for  such  an  integration.  We  recommend  pursuing  this 
integration  with  an  appropriate  engineering  scale  propagation  modeling  system  such  as  WaveTrain  by 
MZA  Associates. 

•  Further  Application  to  Large-Scale  Problems:  We  anticipate  applying  this  new  modeling  capability  to 
aircraft  and  systems  of  interest  to  the  Air  Force  and  other  agencies. 
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